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ABSTRACT 

Asymptotic-induced  methods  are  presented  for  the  numerical  solution  of  hyperbolic  con¬ 
servation  laws  with  or  without  viscosity.  The  methods  consist  of  multiple  stages.  The  first 
stage  is  to  obtain  a  first  approximation  by  using  a  first-order  method,  such  as  the  Godunov 
scheme.  Subsequent  stages  of  the  method  involve  solving  internal-layer  problems  identified 
by  using  techniques  derived  via  asymptotics.  Finally,  a  residual  correction  increases  the 
accuracy  of  the  scheme.  The  method  is  derived  and  justified  with  singular  perturbation 
techniques. 
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1  Introduction 


The  combination  of  asymptotic  and  numerical  analyses  provides  improved  accuracy  for 
multiple-scales  problems*.  Many  physical  problems  have  multiple  scales;  a  typical  situa¬ 
tion  occurs  when  physics  on  the  fastest  scale  induces  narrow  regions  where  the  variation  in 
the  solution  is  large.  Such  regions  are  called  boundary  layers  or  transition  layers,  depending 
on  whether  they  are  near  a  boundary  or  inside  the  interior  of  the  domain.  Examples  of 
such  situations  are  laminar  flow  of  a  slightly  viscous  fluid  or  combustion  with  high  activa¬ 
tion  energy.  Classical  schemes  applied  to  these  types  of  situation  generally  fail  to  correctly 
describe  the  behavior  inside  the  layers.  This  difficulty  is  overcome  by  developing  numerical 
methods  based  on  the  asymptotic  analysis  of  the  multiple-scales  problems.  We  demonstrate 
how  asymptotic  analysis  and  numerical  analysis  can  interact  to  achieve  the  goal  of  a  highly 
accurate  solution  method. 

The  numerical  techniques  developed  in  this  paper  are  designed  to  solve  two  problems: 
hyperbolic  conservation  problems  with  or  without  viscosity.  Also,  we  treat  three  different 
types  of  singularity  arising  in  inviscid  and  viscous  conservation  laws,  using  domain  decom¬ 
position:  specifically,  we  solve  problems  that  have  shocks,  weak  singularities,  and  interaction 
of  singularities.  We  use  an  asymptotic  t  nalysis  to  identify  the  appropriate  scalings  as  well 
as  modified  problems  in  the  transition  layers.  The  domain  decomposition  of  our  numerical 
method  is  then  based  on  the  criterion  given  by  the  asymptotic  analysis.  The  numerical  algo¬ 
rithms  involve  several  stages  to  solve  each  modified  problem  in  the  layers.  Finally,  we  treat 
weak  singularities  with  a  residual  correction.  The  residual  correction  also  identifies  incorrect 
asymptotic  preconditioning  and  controls  the  numerical  accuracy  of  the  solution. 

A  number  of  recent  papers  present  new  numerical  methods  that  are  induced  by  am  asymp¬ 
totic  anailysis.  For  example,  basis  functions  motivated  by  asymptotics  of  boundary  layers 
are  developed  in  [3].  Many  of  the  basic  ideas  relating  to  asymptotic  analysis  and  numerical 
methods  that  use  domain  decomposition  are  found  in  [2].  These  ideas  were  incorporated  into 
a  parallel  numerical  method  in  [15].  Specific  applications  to  problems  involving  conservation 
laws  have  been  developed  in  [1]. 

On  the  other  hand,  some  recent  papers  used  the  singular  perturbation  approach  to  in¬ 
vestigate  the  effect  of  viscosity  on  conservation  laws.  Sharp  estimates  can  be  found,  for 
example,  in  [10]  and  its  references.  The  dynamics  of  internal  layers  are  studied  in  [6].  Also, 
uniform  approximation  based  on  the  matching  asymptotic  technique  are  given  in  [8]  and  [7]. 

This  paper  is  a  further  stage  toward  the  synthesis  of  asymptotic  analysis  and  numerical 
methods  applied  to  systems  of  viscous  conservation  laws  (seen  as  a  singular  perturbation 
problem).  We  begin  with  an  asymptotic  analysis  of  the  problem  in  Section  2.  This  anal¬ 
ysis  is  expanded  and  combined  with  numerical  techniques  in  Section  3.  We  end  with  a 
demonstration  of  the  techniques  on  computational  fluid  dynamics  problems  in  Section  4. 

2  Asymptotic  Analysis 

In  this  section  we  present  the  problem  and  briefly  review  the  asymptotic  analysis 
in  the  method.  The  result  is  a  uniform  approximation  based  on  the  asymptotic  technique  of 
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matching  [5]. 

Asymptotic  analysis  that  is  strongly  coupled  with  the  numerical  analysis  is  delayed  until 
later  in  the  paper. 


2.1  Setting  of  the  Problem  and  Regular  Expansion 

Consider  the  Cauchy  problem 

f  +  £fW)  =  «£  (JWgf)  tol 

i 

t/(z,0)  =  V’(x)  for  x  G  1R. 


(i) 


Here  the  solution  U  :  fi  — »  1R"  is  a  vector- valued  function  ,  the  domain  is  fi  =  IRx]0,T[, 
and  e  <<  1  is  a  small  parameter. 

We  assume  that  V  is  piecewise  smooth.  We  also  assume  that  F  and  P  are  smooth 
functions  of  U.  We  suppose  that  P  is  a  suitable  viscosity  matrix  [4]  for  the  shocks  of  the 
following  associated  inviscid  problem: 


’8  +  W  =  0  for  (*,0  6  0 

U{ i,0)  =  V(x)  for  x  £  IR. 


(2) 


That  is,  a  shock  wave  solution  to  (2)  can  be  obtained  as  a  limit  of  progressive  wave  solutions 
of  (1).  Problem  (1)  is  a  parabolic-hyperbolic  singular  perturbation  problem  driven  by  (2). 
One  easily  obtains  the  regular  expansion 

tcter  =  U°  +  tUl  +  e3t/3  +  . . .  (3) 


that  is  a  priori  valid  outside  the  neighborhood  of  the  singularities  of  the  solution  to  (2).  We 
substitute  U™Ut  in  the  differential  equation  of  (1)  and  use  identification  in  e  to  obtain  that 
U°  must  be  a  solution  of  (2).  We  also  find  that  U1  must  be  a  solution  of  the  following  linear 
hyperbolic  problem: 


es-  +  &(DF{V°)U')  =  ii(P(V°)% £)  for 

f/1(x,0)  =  0  for  x  G  IR, 


(4) 


where  DF  denotes  the  Jacobian  matrix  of  F.  In  general,  the  equation  for  the  term  U*  with 
coefficient  e*  is  the  solution  to  the  analogous  differential  equation 

^  ^  (• DF(U °)IP)  =  RHS(Uj,j  <  i), 

with  a  right-hand  side  that  depends  on  the  previous  terms.  The  inviscid  problem  (2)  has 
many  weak  solutions;  we  shall  uniquely  define  U°  as  the  analysis  progresses. 

It  is  most  common  to  complete  t/^uker  with  a  f/^”er  for  each  singular  region  to  obtain  a 
uniformly  valid  approximation  to  a  viscous  problem.  However,  we  shall  develop  in  detail  only 
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the  asymptotic  analysis  that  is  required  for  our  numerical  methods.  Thus  we  restrict  our 
asymptotic  analysis  to  complete  U^a  with  only  when  the  singularity  is  a  shock  layer. 
The  solution  is  completed  for  weak  singularities  with  the  residual  correction  techniques  of 
Section  3.4.  This  completion  also  results  in  a  uniformly  valid  approximation.  Furthermore, 
the  approximation  is  not  difficult  to  obtain  numerically.  It  is  easy  to  see  that  interaction  cf 
singularities  requires  a  stretching  in  both  the  space  and  time  variables:  this  leads  to  some 
parabolic  layers.  We  shall  use  a  numerical  treatment  of  the  interaction  of  singularities  that 
does  not  require  a  complete  analysis  of  these  parabolic  layers. 


2.2  Transition  Layer  That  Corresponds  to  Shock  Waves 

We  assume  that  solutions  U°  of  (2)  are  smooth  except  on  piecewise  regular  curves  S*(t).  It 
is  assumed  for  t  6  [to  h]  that  the  S*  are  isolated  from  each  other.  Without  loss  of  generality 
we  may  drop  the  subscript  k.  S  is  assumed  to  be  smooth  for  t  €  [£o,fi]-  In  particular,  we 
assume  no  focusing  of  characteristics.  To  make  this  more  precise,  we  suppose  that  there  is 
an  interval  of  time  [to,  ti]  such  that  for  each  t  £  [to,fi]  and  for  each  S  the  following  limits 
exist: 

U?=  lim  U°(x,t),  U°=  lim  U°(x,  t), 

x—S-(0  X— S+(0 

Let  us  define  the  change  of  variable:  x  =  x  —  S'(t)  t  and  r  =  t.  We  denote  U(x,t )  = 
U(x,t).  We  further  assume  that 


dU° 


dU° 


~  &  dr  ’  ^  dr  ’ 


t/‘,  =  lim 
3’  *—o 


ftu3 

Ujr  -=  lim  ,  for  t,  j  >  0, 

4-.0+  C7Z* 


dx*  ’ 

u?,  =  0{  1),  Ul  =  0(l). 

The  shock  layer  profile  will  not  have  rapid  variation,  so  it  is  appropriate  to  scale  and 
translate  only  the  spatial  variable  (and  not  the  temporal  variable).  Such  a  transformation 
is  defined  by 

£  =  - — and  r  =  t,  (5) 


where  we  denote  U((,t)  =  U(x,t).  Under  this  transformation  the  differential  equation  of 
problem  (1)  becomes 


dU_ 

dr 


(F(U)  -  S(t)U) 


(6) 


This  suggests  an  inner  expansion  of  the  form 

=  t/°  +  eU1  +  e3U2  +  . . . ,  (7) 

where  all  of  the  terms  ft'  are  functions  of  (  and  r.  Using  this  expansion  in  (6)  and  imposing 
the  matching  relations  [5,  8]  with  the  outer  expansion  U™itT,  we  derive  the  set  of  ODE 


3 


problems  for  the  O'.  The  equation  for  the  first  term  is 

■  -& +  §i  (no°)  -  mu°) = 

0°  -»  Uf  as  (  -*  -oo 
.  0°  -*  U?  as  (  ->  +oo 


and  for  the  second  term  we  have 

'  ( P(U° )  •  tf1)  +  £  (DF(f/°)  •  i/1  -  S'(t)U')  = 

-  II u1  ~  {U^i  +  t/°,)|U  ->0  as  f  -*  -oo  (9) 

.  HI/1  -  (t/0V(  +  ^.Jllo  -  0  as  i  -  +oo, 
where  ||(/||0  =  mai,=i..n(|ui|)  for  f/  =  (uj)i=i..n.  More  generally,  we  obtain 

'  (W0)  •  #*)  +  it  (DF(U°)  •  -  S'(t) O')  =  RHSi{U’,j  <  t), 

•  as  ^  >  — oo  (10) 

,  ll^-rU^^r^llo-O  as  t  -+  +oo 

where  the  right-hand  side  is  a  nonlinear  function  depending  on  the  03  for  j  <  i.  Identification 
of  the  coefficients  of  the  appropriate  power  of  e  in  the  differential  equation  of  problem  (6) 
determines  RHS,. 

The  temporal  variable  r  can  be  considered  as  a  parameter  in  the  transition  layer  because 
the  above  problems  require  the  solution  of  only  ODEs.  Also,  the  problems  for  O'  are  linear 
for  i  >  0.  Being  able  to  treat  r  as  a  parameter  is  not  surprising,  since  U°t  =  0(1)  and 

U°,T  =  0(1). 

Next  we  discuss  the  existence  of  Oo  and  the  uniqueness  of  Uo  as  well  as  the  uniqueness 
of  the  curve  S.  Using  a  matching  relation  on  the  first  spatial  derivates  on  the  terms  in  the 
expansion  U™a,  one  looks  foi  U°  such  that 


— - >  0  as  (  — ♦  rfcoo. 

oi 

We  integrate  (8)  from  — oo  to  (  to  obtain 

.  d0° 

p(0t  =  «(0-»W),  (11) 

where  H{U )  =  F{0)  —  S'(t)0.  Thus,  the  existence  of  a  solution  to  (8)  implies  the  Rankine 
and  Hugoniot  condition  (RH) 

H{U?)  =  H(U?). 

This  means  that  U°  and  U°  are  critical  points  of  the  dynamical  system  (11). 

In  addition  to  the  RH  condition,  we  must  be  able  to  construct  the  layer;  thus,  we  assume 
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HI. 


There  exists  a  unique  trajectory  for  the 
dynamical  system  (11)  from  Uf  to  U°. 


It  is  interesting  to  compare  HI  with  the  classical  geometric  entropy  condition  (GEC).  For 
scalar  conservation  laws,  HI  is  equivalent  to  the  GEC  [11]  [14];  however,  this  is  not  in  general 
true.  A  number  of  very  interesting  studies  have  considered  which  conditions  one  can  expect 
that  the  GEC  implies  HI.  In  particular,  the  GEC  implies  HI  when  P  is  the  identity  and 
the  problem  consists  of  a  2  x  2  system  with  the  assumption  of  genuine  nonlinearity  [4,  12]. 
The  existence  of  a  viscous  profile  (i.e.,  HI)  is  not  of  the  same  nature  as  the  GEC.  However, 
we  restrict  our  problems  to  cases  where  the  RH  condition  and  HI  are  enough  to  uniquely 
define  U°  and  S  (and  consequently  U™ta). 

Now  we  have  determined  U°  up  to  a  translation  in  the  spatial  variable  In  Section  2.2 
this  translation  will  be  imposed  as  a  solvability  condition  on  the  solutions  to  problem  (9). 


2.3  Construction  of  the  Inner  Expansion  in  the  Shock  Layer 

Here  we  briefly  outline  the  construction  of  the  solution  in  the  shock  layer.  This  formal  theory 
for  systems  of  conservation  laws  is  a  generalization  of  the  construction  in  the  case  of  scalar 
conservation  laws  [8].  For  more  details  of  this  generalization,  see  [9]. 

Let  us  concentrate  on  the  construction  of  the  solution  U 1  of  problem  (9).  The  procedure 
is  first  to  obtain  two  functions  {//(£,  t)  and  U,(£,t)  as  solutions  to  related  problems,  then 
combine  them  such  that  the  resulting  function  is  a  smooth  solution  of  (9).  The  function  17, 1 
satisfies  the  equation  from  problem  (9)  with  the  left  boundary  condition 

Ul  -  (f/oy  +  f/x°()  -  0  as  (  —  -oc, 

whereas  fir*  satisfies  the  equation  from  problem  (9)  with  the  right  boundary  condition 

U1  -  (UK  +  Ul)  ->  0  as  t  -»  +co. 


We  use  the  solutions  to  these  problems  that  are  constructed  using  an  element  dU°/d(,  of 
the  kernel  of  (9).  That  is,  we  construct  Ul  as  the  function 


u'U,t) 


tVU.O  +  B,(t)xSg.  for  £<0 
UXt,t)  +  Br{t)x2g.  for  f>0, 


where  U  x  V  denotes  the  vector  (u,u,)t=1.n.  Bi  and  Br  are  some  smooth  vector  functions 
from  fl  into  IR  that  satisfy  the  continuity  condition 


0?(O,i)  +  Bi(<)  x  -g£-(°,0  =  0?(O ,t)  +  Br(t)  x  —  (0,*)- 

Notice  that  for  a  fixed  value  of  t,  this  relation  determines  the  vector  functions  up  to  a  single 
constant. 
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One  can  prove  [9]  that  U 1  is  a  smooth  function  iff  U°  satisfies  the  relation 

I  {/I  t]  -  d(+r  i)  - 


5f/°i 

=  S'(t)Ul  -DF(U°)-Ul  +  P(U0)—  , 

where  [■]  denotes  the  jump  across  the  shock.  This  is  simply  an  area  relation  that  determines 
the  shift  in  (  for  U°.  Let  us  assume  for  simplicity  that  U  and  its  space  derivatives  vanish  at 
infinity.  One  can  show  that  this  relation  is  a  consequence  of  the  conservation  relation 


satisfied  by  solutions  to  (1)  and  (2).  This  relation  is  also  be  satisfied  by  our  uniform  approx¬ 
imation  of  the  solution  to  (1),  where 

f/M  =  (1  -  H(tu))UT“  +  +  eU1)  +  0(e2). 


Here  =  (i  —  S(t))/ev  is  the 


intermediate  variable  and  H  is  a  smooth  cutoff  function 


if  | y|  >  2 

if  \y\  <  1. 


More  precisely,  we  have 


dUa± 

7-oo  dt 


=  0(e), 


We  currently  have  defined  Ul  up  to  a  kernel  function  of  (9).  Since  we  use  only  the 
first  term  U°  in  the  current  implementation  of  this  method,  we  do  not  discuss  how  to 
determine  this  constant  here.  This  construction  can  be  pursued  at  any  order.  For  example, 
the  solvability  condition  for  problem  (10)  with  i  =  2  will  select  the  only  admissible  U 1 
solution  of  problem  (9)  [8]. 


3  Asymptotic-Induced  Numerical  Scheme 

Now  we  consider  the  numerical  treatment  of  problems  (1)  and  (2)  by  using  the  analysis  in 
the  previous  sections. 


3.1  Localization  of  the  Singularities 

Let  {*<}»=..., -1,0,1,...  be  the  spatial  discretization  of  step  size  Ax,  and  let  be  the 

temporal  discretization  of  step  size  At.  For  convenience  we  restrict  the  discussion  on  local¬ 
ization  of  singularities  to  scalar  conservation  laws.  The  scalar  terms  are  distinguished  from 
the  system  by  using  the  notation  u  and  /  as  the  solution  and  flux  function  in  place  of  U  and 
F,  respectively.  We  suppose  that  u  is  known  at  time  t*.  Starting  from  u(-,t*),  we  apply  a 
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Table  1:  Asymptotic  Order  of  Residual 


Type  of  Zone 

Order  of  Residual 
+  /(u)* 

Local  Coo 

£ 

rdinates 

r 

Regular  zone 

0(e) 

X 

t 

Shock  layer 

0(e~') 

(x  -  S(t))/e 

t 

Weak  singularity 

0(7 7*) 

(x  -  S(i))/«,/2 

t 

Shock  interaction 
with  other  singularities 

OfP1)- 

(x  —  So  -  S\t)/e 

(t  -  to)/* 

Discontinuous  with 
/  locally  linear 

0(1) 

t  -  ta 

Formation  of  shock 

0(e~1/*) 

discretization  that  gives  an  approximation  uK  at  the  discrete  spatial  points  {x,}  at  time  , 
where  K  >  k.  Given  this  approximation,  we  discuss  the  problem  of  localization  of  different 
types  of  singularities  of  u  at  time  t#. 

Localizing  the  singularities  will  use  the  residual  obtained  from  the  left-hand  side  of  the 
PDE  in  (1).  To  this  end,  we  recall  in  Table  1  some  results  of  the  asymptotic  analysis  of  the 
Cauchy  problem  (1);  see  [8]  and  its  references.  Here  S0  and  Si  are  the  constants  from  the 
local  linearization  S(t )  =  So  +  Sit  +  0((t  —  ta )2)  in  the  neighborhood  of  the  singularity. 

The  localization  is  presented  for  the  case  when  uK  is  obtained  by  using  the  first-order 
Godunov  scheme;  however,  this  analysis  requires  only  minor  modification  for  several  first- 
order  discretizations  of  scalar  conservation  laws,  provided  the  discretizations  are  conservative. 
The  solution  obtained  via  the  Godunov  scheme  is  an  approximation  to  within  0(Ax2)  of  the 
following  Cauchy  problem: 


fH  +  £/(«)  =  Ax£  (G  (/,  u,  At,  Ax)  §h)  t6[tk,  tK } 
u(x,tk)  given. 


(12) 


We  fix  the  ratio  At/Ax  so  that  Ax  may  be  treated  as  the  parameter  e  in  the  table. 

In  the  more  general  case  of  a  system  with  some  first-order  discretization,  we  obtain  Wq 
as  the  numerical  approximation  of  the  solution  to  the  hyperbolic  problem.  Let  us  assume 
that  the  function  Wg  is  the  solution  to 

'  f  +  £F(V)  =  Az£(G„(F,V,A*,A!)g)  I  6  [(*,(*] 

< 

.  U(x,tk)  given, 


where  G#  is  for  the  numerical  viscosity  of  the  specific  hyperbolic  scheme  with  Gh  =  0(1). 
The  detection  for  the  numerical  method  is  based  on  obtaining  an  approximation  to  the 
residual  dU/dt  +  dF(U)/dx. 

The  residual  is  of  magnitude  0(Ax_1)  in  either  a  shock  layer  or  in  a  zone  where  a  shock 
interacts  with  some  other  singularities.  In  addition,  an  approximation  of  the  viscous  term 
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■§£jU(-,tK)  localizes  some  of  the  singularities.  For  example,  this  viscous  term  will  be  of  order 
0(Ax~1)  in  a  shock  layer  or  in  a  zone  of  interaction. 

Numerically  it  may  be  difficult  to  determine  when  the  residual  (or  viscosity)  is  of  order 
Ax-1.  This  difficulty  may  usually  be  overcome  by  computing  the  numerical  solution  on 
two  grids  with  different  spatial  spacings  Aii  and  Axj.  Let  R\  and  R2  be  approximations 
to  the  residuals  of  the  solutions  on  these  meshes.  When  the  ratio  Ri/Rj  of  the  residuals  is 
within  some  tolerance  of  Axj/Axj,  then  we  treat  the  region  as  being  inside  a  shock  (or  shock 
layer).  It  should  be  mentioned  that  asymptotic  analysis  and  numerical  experiments  show 
under  general  conditions  that  strong  singularities  are  not  well  localized  by  any  approximation 
of  the  gradient  of  the  solution. 

Shocks  (and  interactions  with  shocks)  are  easily  detected  and  are  treated  in  the  next 
section.  In  Section  3.4,  we  shall  show  how  some  weak  singularities  that  are  much  more 
difficult  to  localize  can  be  absorbed  in  a  residual  correction  process. 


3.2  Numerical  Treatment  of  Shocks 

Let  us  first  recall  the  numerical  treatment  of  a  simple  shock  for  a  scalar  conservation  law  [1], 
Let  fl0  =  [<*,  &]  be  the  zone  that  includes  the  shock  (or  the  shock  layer).  First  we  construct 
an  approximation  w  to  the  solution  of  (2);  next,  to  compute  an  approximation  to  the  shock 
location  5,  we  treat  the  solution  near  the  shock  (or  inside  the  shock  layer). 

The  region  fi0  is  chosen  so  that  u  is  an  a  priori  valid  approximation  to  u  except  possibly 
for  some  portion  of  fi0;  thus,  w  =  u  outside  Qq.  The  approximation  w  is  constructed  inside 
with  an  extrapolation  of  u.  This  extrapolation  is  based  on  values  of  u  outside  Q0.  We 
construct  w  as  the  connection  of  the  extrapolated  values.  For  t  —  tK  fixed,  we  impose  the 
area  relation 


on  the  discrete  problem.  This  is  presented  visually  in  Fig.  1,  where  we  use  second-order 
extrapolation  and  have  placed  S  so  that  the  area  covered  by  the  dark  grey  is  equal  to  the 
area  covered  by  the  light  grey.  As  we  discuss  below,  the  left  and  right  values  of  w  are  used 
for  the  boundary  conditions  for  the  first-order  approximation  of  the  solution  with  a  viscous 
profile. 

Notice  that  the  accuracy  of  the  approximation  w  as  a  solution  to  the  inviscid  problem 
is  limited  only  by  the  particular  numerical  techniques.  Specifically,  the  accuracy  is  limited 
by  the  accuracy  to  which  u  is  computed  outside  fi0  and  by  the  order  of  accuracy  of  the 
extrapolation.  There  is  no  limitation  imposed  by  the  asymptotic  analysis  on  this  aspect  of 
the  method. 

The  treatment  of  the  solution  in  a  neighborhood  of  the  shock  (or  inside  the  shock  layer) 
depends  on  whether  the  goal  is  the  solution  of  the  viscous  problem  (1)  or  the  inviscid  problem 
(2).  If  the  solution  to  inviscid  problem  is  the  goal,  then  we  wish  to  minimize  the  viscosity. 
Viscous  problems  require  the  computation  of  the  profile  of  the  shock.  The  asymptotic 
analysis  of  the  shock  layer  for  a  system  of  conservation  laws  is  similar  to  the  case  of  a 
scalar  conservation  law  (cf.  Section  2). 

For  the  numeric  computations,  the  primary  difference  between  the  treatment  for  a  scalar 
conservation  law  and  a  system  is  that  we  must  verify  that  certain  conditions  are  satisfied  oy 
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Figure  1:  Approximation  of  shock  location 


the  system.  Specifically,  we  a  posteriori  (during  the  numerical  computations)  verify  that 

•  W°  and  Wj°  are  on  a  Rankine-Hugoniot  set  up  to  a  tolerance, 

•  there  is  a  trajectory  from  W°  to  W°,  and 


=  0(1). 


The  first  assumption  is  satisfied  automatically  by  solutions  to  scalar  conservation  laws.  Now 
we  shall  describe  in  more  detail  the  numerical  treatment  of  the  shock  layer  depending  on 
whether  one  solves  an  inviscid  or  a  viscous  problem. 


3.2.1  Minimum  Viscosity  Method 

Suppose  that  for  a  particular  <*,  the  shock  location  5  is  in  [®j,  x^+1  [,  where  [x,,  x,+i  [C  fio-  In 
the  most  general  case,  5  will  not  be  a  grid  point,  and  we  modify  W  to  maintain  conservation. 
This  step  results  in  modifying  either  the  value  of  Wf  or  the  value  of  Wf+1  (see  [1]  for  details). 
The  modification  introduces  the  minimum  of  viscosity  needed  to  use  the  Godunov  method 
throughout  the  domain.  This  modification  of  Wf  or  is  the  final  stage  of  the  shock  layer 
treatment  when  the  goal  is  the  solution  to  the  inviscid  problem. 

3.2.2  Viscous  Problems 

To  solve  problem  (1),  we  implement  the  inner  layer  computation  as  in  Section  2.2.  The 
matching  relations  used  to  obtain  a  uniform  approximation  are  based  on  the  numerical 
approximation  of  the  outer  expansion.  Thus,  we  solve  (8)  except  that  W°  and  W,0  are  used 
in  place  of  U°  and  U°,  respectively. 
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To  obtain  the  higher-order  corrections  U'  for  i  >  0,  we  would  solve  the  problems  (9)  and 
(10)  with  modified  boundary  conditions;  however,  an  easy  approximation  can  be  obtained 
by  computing  a  poor  approximation  to  Ul  such  that 


Ul 


W  —  Uf  for  x  <  Xj 
W  —  U®  for  i  >  Xj. 


The  numerical  treatment  of  this  problem  involves  using  an  OPp  solver  to  obtain  the 
solution  to  the  system  of  ODEs  (8).  Notice  that  except  when  the  details  of  the  viscous 
profiles  are  desired,  we  need  only  obtain  the  values  of  the  inner  expansion  at  the  coarse  grid 
points  on  which  we  know  W.  Even  when  a  higher-order  numerical  method  is  used,  there  are 
a  large  number  of  steps  in  the  integration  of  the  ODE  system  for  the  inner  region.  The  cost 
of  this  integration  may  be  reduced  by  using  an  explicit  formula  of  approximation  given  by 
the  asymptotic  analysis  in  the  neighborhood  of  the  critical  points  (see  [1]). 


3.3  Interaction  of  Singularities 

In  regions  where  solutions  to  problem  (1)  contain  shocks  that  interact  with  other  singularities, 
we  use  a  brute-force  approach.  The  local  scaling  in  both  the  spatial  and  the  temporal 
variables  is 

>  x  —  So  i  —  t0 

*■  =  — —•  r  =  ~- 

Under  this  transformation  the  PDE  that  governs  the  solution  becomes 


W  ,  9  pan  -  d  fP(ind0 

+  TzF{u)  ~  dt  {'  ( 


(13) 


where  r)  =  U[x,t).  This  is  the  equation  that  is  solved  in  the  regions  with  interactions. 
There  is  no  expansion  of  the  solution  used  here.  The  transformation  a  priori  resolves  all  of 
the  physics.  This  situation  is  reflected  by  all  of  the  terms  in  (13)  being  order  unity. 

The  initial  condition  for  this  problem  is  derived  by  imposing  C°  continuity.  This  was 
sufficient  for  the  problems  considered  in  this  paper.  However,  the  computation  of  the  residual 
is  a  check  of  the  numerical  accuracy  of  the  composite  scheme.  A  residual  correction  similar 
to  the  one  in  the  next  section  could  be  applied  to  improve  the  result.  Interface  boundary 
conditions  is  a  topic  of  further  study. 


3.4  Residual  Correction 

Second-order  accuracy  in  both  space  and  time  is  achieved  through  a  first-order  numerical 
approximation  followed  by  a  residual  correction.  This  accuracy  is  obtained  in  regions  where 
the  regular  expansion  (3)  is  valid.  We  also  discuss  the  residual  correction  in  the  presence 
of  weak  singularities,  where  the  accuracy  is  increased  (but  not  as  much  as  for  the  smooth 
regions).  This  discussion  begins  with  the  asymptotic  analysis,  is  followed  by  the  choice  of 
the  particular  discretization,  and  ends  with  the  treatment  of  weak  singularities. 
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3.4.1  Asymptotic  Analysis 


The  asymptotic  analysis  is  based  on  the  regular  expansion  presented  in  Section  2.1.  Thus, 
this  analysis  is  restricted  to  regions  where  the  solution  is  smooth  and  the  regular  expansion 
is  valid.  For  convenience,  our  analysis  assumes  e  and  A t  are  both  O(Ax).  The  dependency 
of  e  on  Ax  is  related  to  the  order  of  the  regular  expansion  and  the  hyperbolic  scheme  used. 

The  Godunov  method  applied  to  (2)  gives  us  Wq  as  a  first  approximation  to  the  solution 
to  the  parabolic  problem  (1);  thus  U  =  Wq  +  0(Af,  Ax).  A  single  correction  We  will  be 
computed  so  that 

U  =  WG  +  hWc  +  h*W3  +  0(h3),  (14) 

where  h  is  a  small  parameter  of  magnitude  O(At)  =  0( Ai).  The  equation  for  We  is  deter¬ 
mined  with  an  analysis  similar  to  that  of  the  regular  expansion. 

The  viscous  terms  now  arise  from  both  the  viscosity  in  (1)  and  from  the  numerical  error 
in  WG,  resulting  in  the  linear  problem 


’  +  &(DF(U)WC)  =  R(Wc) 

We(x,  0)  =  0 


(15) 


for  Wc.  Here, 

-  4  -  (tt  +  lF(Wc))  ■  (16) 

where  v  =  e  or  v  —  0,  depending  on  whether  we  are  solving  the  viscous  problem  or  the 
inviscid  problem.  Notice  that  the  terms  of  higher  order  would  be  included  in  the  next  term 
of  the  regular  expansion  of  U .  The  right-hand  side  in  the  equation  of  (16)  is  the  residual 
from  using  Wc  as  an  approximation  to  the  solution  of  (1).  To  obtain  a  second-order  method, 
we  start  with  the  first-order  approximation  WG ,  then  compute  the  right-hand  side  of  (16)  to 
second-order  accuracy,  and  finish  with  a  first-order  scheme  for  (15). 

The  asymptotic  techniques  are  similar  in  the  more  general  case  when  e  is  not  of  size 
0( Ax),  and/or  when  a  different  order  method  is  used  to  compute  WG-  Consider  the  case 
when  the  numerical  method  is  0(Axp,  At9).  In  this  case  we  choose  the  parameter  h  in  the 
expansion  (14)  to  have  magnitude  equal  to  the  largest  of  Axp,  Atq,  and  v.  The  problems 
for  the  terms  in  the  regular  expansion  are  obtained  by  identification  in  h  after  substitution 
of  the  expansion  into 

f + iF{u)  -  4  (p(c,>§) + A7FP'<-U) + At,p'{v)' 

where  v  =  0  or  e,  and  P\  and  P3  are  terms  of  0(1).  We  do  not  require  specific  forms  for  Pi  and 
P3  for  the  numerical  method  since  we  may  use  a  residual  in  place  of  AxpP\{U )  +  AtqP3{U). 


3.4.2  Discretization  of  Correction 

Here,  we  choose  discretizations  for  the  numerical  approximation  of  the  residual  and  the 
correction  scheme.  Notice  that  these  discretizations  as  well  as  the  choice  of  the  Godunov 
method  for  the  first  approximation  are  only  examples  of  possible  treatments. 
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Centered  differences  are  used  throughout  this  section;  hence,  the  resulting  finite-difference 
scheme  will  be  a  conservative  discretization. 

Residual.  Here  we  use  the  discretization  based  on  centered  differences  and  averaging. 
The  temporal  derivative  is  approximated  by  dWa{x%,  as  where  we  use 

the  notation  6  to  be  the  centered  difference  operator 

c  9(V  +  At?)  -  g( t?  -  At?) 

Sl"sM  = - 2 5? - • 

The  formula  for  the  complete  residual  is 

^+1/J  =  l  {-  [w£1/J  +  \  (6>*F*k+1  +  (17) 

[*,(#*."&)  +  4(^+14^l)J }  • 

Here  we  use  the  notation  F*  =  F{Wqx)  and  the  average  Pf+i/?  =  (/>(W'g.J  +  7,(^G,i+i))/2- 
Taylor  series  analysis  can  be  used  to  verify  that  this  is  an  0( Ax2,  At2)  approximation  to  the 
residual  at  the  point  (x,  t)  =  (xl,t/e+ 1/2). 

Correction  Scheme.  The  implicit  scheme  based  on  a  centered  spatial  derivative  and  a 
centered  temporal  derivative  is  used  for  the  discretization  of  Wc  to  obtain 

W*+1/2  +  8lx  ( DF(W SJ1)  •  ^+1)  =  ^*+1/2.  (18) 


This  is  a  first-order  accurate  method. 

The  accuracy  of  the  discretization  for  the  residual  is  balanced  with  the  accuracy  of  the  dis¬ 
cretization  for  the  correction  term.  We  compute  Wq  to  first  order.  The  residual  is  computed 
to  second  order;  however,  it  is  scaled  by  1  jh  for  the  computation  of  the  correction.  This 
results  in  a  first-order  approximation  to  the  source  term  for  the  correction.  The  correction 
discretization  is  first  order;  hence  the  overall  method  is  second  order. 


3.4.3  Weak  Singularities 

A  modification  of  the  above  analysis  can  be  used  to  develop  the  residual  correction  in  the 
presence  of  weak  singularities  because  the  residual  is  still  small  in  the  neighborhood  of  a 
weak  singularity  (cf.  Table  1).  The  numerical  algorithm  is  the  same  for  this  case;  however, 
the  analysis  (and  resulting  accuracy)  are  different.  In  place  of  (14)  we  have  the  regular 
asymptotic  expansion 

U  =  WG  +  W1  +  ..., 

wVere  Wg  is  the  approximation  of  U  obtained  from  the  numerical  method.  The  function 
Wg  is  an  approximation  to  the  solution  of 


dU  d  d 

dt  +  dxF ^  ~vdx 


where  v  =  e  or  v  =  0.  However,  Wg  may  be  treated  as  the  exact  solution  of  the  modified 
equation 

-sf  +  -F(Wo)  =  UlHS,  (19) 
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where  h  =  o(l)  is  defined  such  that  RHS  =  0,(  1)  is  true  a  priori.  So  we  rescale  Wx  as: 
W\  =  hWc.  As  in  Section  3.4.1,  we  obtain  the  equation  for  Wc  as 


dWc 

dt 


+  ftiDF(W„).W')  =  -RHS+l‘L 


This  is  the  equation  for  the  correction  term  in  the  presence  of  weak  singularities. 

The  numerical  discretizations  used  for  the  case  of  the  weak  singularity  can  be  the  same  as 
those  presented  in  the  last  section.  For  example,  when  the  singularity  is  of  order  h  =  0(y/e), 
then  the  method  with  the  same  discretizations  will  have  accuracy  0(h2)  =  0(e).  Both  forms 
of  the  analysis  in  this  section  are  valid  with  the  numerical  discretizations  of  the  previous 
two  sections.  The  numerical  method  after  the  residual  correction  will  be  more  accurate 
than  using  the  first  approximation  Wg  alone.  In  addition,  it  is  feasible  to  use  a  higher-order 
method  for  the  discretizations  of  the  residual  and  a  higher-order  expansion  for  the  correction 
in  regions  of  weak  singularities. 


3.5  Numerical  Algorithm 

The  numerical  and  asymptotic  analysis  are  combined  and  summarized  in  Algorithm  1  below. 
The  outer  region  is  discretized  by  using  the  coarsest  grid.  Numerical  values  of  Wc  and  Wc 
are  determined  at  coarse  grid  points.  Wc ?  is  computed  over  the  whole  grid;  however,  the 
coarse  grid  values  located  inside  refined  regions  are  used  only  in  the  conservation  relation. 
The  refined  grid  values  are  injected  into  these  coarse  grid  values  at  the  end  of  the  time  step. 
Linear  interpolation  is  used  to  obtain  values  of  the  coarse  grid  solution  between  points  when 
necessary. 

The  numerical  results  presented  here  involve  only  a  two-level  refinement.  This  is  because 
the  current  implementation  is  designed  for  problems  in  which  shocks  (or  shock  layers)  and 
their  interactions  are  the  only  strong  singularities.  We  recall  that  weak  singularities  are 
treated  by  the  residual  correction  scheme. 

The  numerical  solution  in  the  shock  layer  is  computed  in  the  transformed  coordinates 
(5)  when  the  shock  profile  is  not  rapidly  varying.  When  interactions  exist,  the  coordinates 
use  the  same  (  with  the  temporal  scaling  r  =  f/e.  Although  A£  is  the  same  magnitude  as 
Ax,  local  spatial  scaling  is  used.  This  means  that  a  region  of  size  Ax  in  the  coarse  grid 
corresponds  to  many  points  in  the  £-grid.  For  this  reason,  we  refer  to  the  £-grid  as  the 
refined  grid. 

The  method  is  adaptive.  The  refined  grid  is  allowed  to  change  shape  and  location  as  the 
solution  evolves.  The  refined  region  contains  an  overlap  region  where  both  the  inner  and 
outer  solutions  are  valid. 


4  Experiments 

In  this  section  we  demonstrate  the  techniques  developed  within  this  paper.  The  experiment  in 
Section  4.1  demonstrates  the  residual  correction.  The  experiments  in  Section  4.2  demonstrate 
the  method  on  the  isentropic  gas  dynamic  equations. 
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For  k  =  1 . 

I.  March  from  f*  to  tk+i  on  two  coarse  meshes  with  spatial  discretizations  Axj  ^ 
Ax2. 

II.  Detection. 

A.  On  each  coarse  mesh  compute  the  residual  used  in  Table  1.  (Alternatively, 
compute  the  viscosity.) 

B.  Mark  regions  that  should  be  refined. 

III.  Compute  Wc  with  the  residuals  of  Step  II. A.  except  in  marked  regions. 

A.  Modify  the  shape  of  the  refined  region. 

IV.  March  the  marked  regions  from  t*  to  t*+i  according  to  the  type  of  singularity. 

A.  For  a  simple  shock  of  Section  2.3,  compute  the  solution  to  the  ODEs  as 
outlined  in  that  section. 

B.  For  an  interaction  of  singularities  of  Section  3.3: 

1.  Form  the  initial  condition  in  newly  refined  regions. 

2.  Determine  boundary  conditions  from  linear  interpolation. 

3.  March  inner  solution  Af/(eAr)  time  steps. 

V.  Inject  the  values  from  the  internal  layer  to  obtain  the  composite  uniform  approx¬ 
imation. 

Algorithm  1  Numerical  algorithm 
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4.1  Solving  Burgers’  Equation 

The  residual  correction  technique  is  demonstrated  by  solving  the  inviscid  Burgers’  equation 

du  du 

Tt+uTz  =  ° 

on  the  domain  [—0.2, 1]  with  a  wedge  initial  condition.  The  first  approximation  is  obtained 
with  a  first-order  Godunov  method  using  discretization  parameters  Ax  =  .025  and  At  = 
.0125.  The  solution  at  time  t  =  .5  is  presented  in  Fig.  2.  The  increase  in  accuracy  of 


X 

Figure  2:  Exact  solution  to  Burgers’  equation 


the  residual  correction  over  using  just  the  Godunov  is  shown  in  Figs.  3-4.  The  minimum 
viscosity  method  is  used  for  the  numerical  solution  in  a  neighborhood  of  the  shock.  A  cutoff 
function  is  used  for  the  residual  correction.  This  causes  the  correction  to  be  applied  only 
outside  a  neighborhood  of  the  shock.  The  mass  that  might  be  lost  or  gained  in  this  procedure 
is  included  in  the  computations  for  the  minimum  viscosity  method.  The  large  jump  at  the 
shock  is  a  reflection  of  the  error  between  the  Li  projection  of  the  numerical  solution  with 
the  exact  solution. 

4.2  Solving  the  Isentropic  Gas  Dynamic  Equations 

In  this  section  the  method  is  demonstrated  on  the  system 

du  dv 
dt  dx 

dt  dx  \u"|7  dx  \dx  J 

Here  u  is  the  inverse  of  the  density  and  v  is  the  velocity.  These  equations  are  obtained 
from  the  conservation  of  mass  and  momentum  in  Lagrangian  coordinates  assuming  that  u 
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is  equal  to  the  pressure  raised  to  the  —  l/7th  power  (the  perfect  gas  law)  with  7  =  2.2.  The 
numerical  solution  in  the  outer  region  is  computed  with  a  first-order  Godunov  method  with 
CFL  number  At/ Ax  =  .2.  The  discretization  on  the  scaled  coordinates  inside  the  shock  layer 
is  based  on  A£  =  .1.  We  use  a  third-order  Runge  Kutta  method  to  compute  the  viscous 
profile  for  a  shock  layer.  To  compute  a  parabolic  layer,  we  use  an  explicit  first-order  ENO 
scheme  (similar  to  the  Godunov  scheme)  with  the  CFL  condition  At/A(  <  .2  and  stability 
condition  At/A( 2  <  .25.  These  values  are  within  the  limits  imposed  for  the  stability  of  the 
finite  difference  method. 


4.2.1  Simple  Shock 


The  first  experiment  is  a  simple  shock.  In  Fig.  5  we  show  the  computed  solution  using 
various  methods.  The  initial  condition  is 


where 


u(i,0)  = 


t 


v(x, 0)  - 


UR  =  2.50,  UL  = 


and 


VL  =  VR  +  \ 


Ui,  for  x  <  0 

UR,  for  x  >  0 

(20) 

VL,  for  x  <  0 

Vr  ,  for  x  >  0 

(21) 

.800,  VR  =  .600, 

riae  value  for  Vi  is  chosen  using  the  Rankine-Hugoniot  condition 


Vl-Vr  _  l/UZ'l/VZ 
Ur  -Ul  Vr~  Vl 


The  computations  are  run  with  Ai  =  .02.  The  Godunov  and  minimum  viscosity  metnods 
approximate  the  solution  to  the  hyperbolic  problem,  whereas  the  ODE  layer  (solving  problem 
(8)  inside  the  shock  layer)  and  the  parabolic-layer  methods  are  used  to  compute  the  numerical 
approximation  to  the  viscous  problem  with  e  =  .01. 

The  minimum  viscosity  method  introduces  significantly  less  viscosity  than  the  Godunov 
method.  As  expected,  the  parabolic- layer  and  the  ODE-layer  methods  produce  numerical 
solutions  that  are  very  close  to  each  other. 


4.2.2  Shock-Rarefaction  Interaction 


The  method  was  also  used  to  solve  a  Riemann  problem  with  a  right-traveling  shock  and 
a  left- propagating  rarefaction.  Thus,  for  the  small  time,  we  must  solve  an  initial  layer 
with  a  shock-rarefaction  interaction.  We  can  compute  the  exact  solution  of  the  inviscid 
problem.  First,  an  analytic  self-similar  solution  (a  rarefaction  emanating  from  the  origin)  to 
the  inviscid  isentropic  gas  dynamic  equations  is  given  by 


u(x,t)  =  71/'('r+1) 


(22) 
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Plot  of  U  showing  various  methods,  t«.4 

Figure  5:  Numerical  solutions  to  the  isentropic  system,  t  =  .4 

v(x, t)  =  ^  -b  const.  (23) 

An  initial  condition  with  a  shock  and  rarefaction  emanating  from  the  origin  is  constructed 
by  connecting  left  values  (f4,  14)  to  middle  values  ( U0 ,  VD)  with  a  rarefaction.  The  middle 
values  are  connected  to  the  right  values  ( UT ,  VT)  with  a  shock.  In  our  experiment  the  initial 
condition  is  given  by  (20), (21)  where 

UL  =  1.4709,  UR  =  2.5000,  14  =  1-0388,  VR  =  0.8050. 

The  middle  value  of  the  solution  between  the  shock  and  rarefaction  is  (U0,V0)  =  (1.973, 
1.356).  We  expect  the  viscous  perturbation  to  have  little  or  no  effect  on  the  speed  at  which 
shocks  and  rarefactions  travel;  thus,  we  shall  compare  the  viscous  solutions  to  the  exact 
solution  of  the  inviscid  problem  given  above. 

In  Fig.  6  the  exact  solution  to  the  inviscid  problem  is  compared  with  the  computed 
solution,  where  Ax  =  .02  and  e  =  .01.  This  plot  shows  that  the  shock  speed  and  rarefaction 
propagation  are  nearly  the  same  for  both  solutions.  Greater  detail  is  shown  for  the  shock 
layer  in  Fig.  7.  Here  we  use  6X  —  2e.  This  plot  demonstrates  that  the  method  has  the  correct 
shock  speed  and  that  the  solutions  are  converging  to  the  inviscid  solution  as  e  j  0.  Detail  of 
the  rarefaction  for  the  same  runs  is  shown  in  Fig.  8.  Here  the  offset  to  the  location  of  the 
rarefaction  is  Ax/2  and  can  be  attributed  to  the  initial  condition  for  the  parabolic  layer. 

This  is  a  domain  decomposition  method.  The  internal-layer  domain  is  detected  by  com¬ 
paring  the  second  derivate  to  tolerance  of  conste.  The  jump  in  internal-layer  subdomain 
boundary  as  depicted  in  Fig.  9  is  caused  by  the  steady  smoothing  of  the  rarefaction  until 
it  is  no  longer  detected  by  the  residual  or  second  derivative  test.  Also  for  t  >  0.2  a  simple 
shock  is  detected  and  can  be  solved  by  the  ODE  layer  method. 
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lnviscid  U  Comparison  With  Numeric 


Figure  6:  Viscous  numeric  compared  with  inviscid  exact,  t  =  .2 


Convergence  to  Invisicid  Rarefaction 
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